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Abstract 

This paper formulates the general methodology for 
estimating the bias error distribution of a device in a 
measuring domain from less accurate measurements when 
a minimal number of standard values (typically two 
values) arc available. A new perspective is that the bias 
error distribution can be found as a solution of an intrinsic 
functional equation in a domain. Based on this theory, the 
scaling- and translation-based methods for determining the 
bias error distribution are developed. These methods are 
virtually applicable to any device as long as the bias error 
distribution of the device can be sufficiently described by a 
power series (a polynomial) or a Fourier series in a 
domain. These methods have been validated through 
computational simulations and laboratory calibration 
experiments for a number of different devices. 

1. Introduction 

The measurement of an instrument always has error 
defined as the difference between the measured value and 
the true value. The total error is the sum of the bias error 
and the random error. Most books and articles on 
uncertainty analysis have studied in detail statistical 
estimates of the random error [1,2]. However, the bias 
error, which is the fixed, systematic component of the total 
error, is not sufficiently discussed because it is often 
assumed that all bias errors have been eliminated by 
calibration. Indeed, the bias error can be determined by 
comparison with a standard having accuracy much better 
than the device being tested. The standard device is 
ultimately traceable to a national or international standard 
laboratory. If the standard values over the whole range of 
measurements are known, the determination of the bias 
error distribution of an instrument is extremely trivial. 
Unfortunately, it is not always possible to have a standard 
device available for calibration. Also, the standard device 
itself has limited accuracy. Finley [3] proposed a unique 
idea for extracting the absolute bias error of an angular 
measurement device by comparing it with another 
instrument of comparable quality. The differences in 
readings from the two devices are obtained with two 


t Research Sciemisl. Model Systems Branch, Member AIAA 
Copyright €> 2001 by the American Institute of Aeronautics and 
Astronautics, Inc. No copyright is asserted in the United States under 
Title 17, U.S. Code. The U.S. Government has a royalty-free license to 
exercise all rights under the copyright claimed herein for Governmental 
purposes. All other rights are reserved by the copyright owner. 


different initial angles of the devices. Next, a Fourier 
analysis of ihc two sets of the differences recovered the 
bias errors for both devices. Finley’s method was further 
discussed by Snow [4] from the standpoint of the Fourier 
transform. Finley’s work show's that in angular 
measurements the absolute bias error distributions can be 
obtained by using two less accurate devices. Naturally, a 
legitimate question is whether there is a general method for 
estimating the bias error distribution of a device. This 
question will be answered in this paper. 

This paper formulates the general methodology for 
estimating the absolute bias error distribution of a device 
in a measuring domain when only a few standard values 
(typically two values) arc known. In the proposed 
approach, the bias error distribution is sought as a solution 
of an intrinsic functional equation in a measuring domain. 
The scaling-based method and translation-based method 
are developed, in which the analytical solutions to the 
functional equations for the bias error distributions arc, 
respectively, expressed as a power series and a Fourier 
series in a given domain. The scaling-based method is 
effective for a large class of the bias error distributions that 
can be sufficiently described by a power series or a 
polynomial. The approximate scaling-based method for 
practical implementation is developed, requiring two 
standard values to determine the complete bias error 
distribution. An empirical rule for selecting an appropriate 
order of a polynomial is suggested to achieve the good 
accuracy of calculating the bias error distribution. Effects 
of the random measurement perturbations on calculation of 
the bias error distribution are studied through Monte Carlo 
simulations. In contrast to the scaling-based method, the 
translation-based method has relatively limited 
applications, but it is particularly useful for angular 
measurements since the bias error distribution can be 
naturally expressed as a Fourier series. The scaling- and 
translation-based methods for simultaneously estimating 
the bias error distributions of two devices arc also 
described. Computational simulations and laboratory 
experiments for calibrating a number of different devices 
are conducted to validate the proposed methodology to 
recover the bias error distribution. 

2. Bias Error and Intrinsic Functional Equation 

Let ni( x ) be the measurement value of a ‘true’ 
physical quantity x by a device having a deterministic bias 
error £(x) then, 
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m(x) = x + £(x). ( 1 ) 

The symbol m(x) can be interpreted as a measurement 
operator of the variable x. In order to obtain an additional 
independent equation for £( x ) , a quantity ax+ J3 that is 
a linear transformation of x is measured by the same 
device, where a is a scaling constant and P is a 
translation constant. The measurement value of a x + P is 
expressed as 

m( ax + P ) = ax + P + £( ax + fi). ( 2 ) 

Eliminating x in Eqs. (1) and (2), one obtains a functional 
equation for e( x) 

S(x) = £(x)-a~ I £(ax + ft ), (3) 

where the difference S(x) is a known function that can be 
measured by the device for given a and P , i.e., 

S( x ) = m( x)-a~ l ( m( ax + P)-P). (4) 

Eq. (3) is an intrinsic functional equation governing the 
bias error distribution £( x). Because no assumption has 
been made, this functional equation for the bias error 
distribution is virtually applicable to any measurement 
instrument. Although more complicated functional 
equations can be similarly constructed as a model of the 
bias error, Eq. (3) enjoys its simplicity without loss of 
generality. Now, the fundamental problem is to find a 
solution to Eq. (3) for £(x) in a given domain 

D = [x } ,x 2 ] and a class of admissible functions 

(including domains and ranges). The discussions on 
functional equations from a mathematical perspective can 
be found in references 5 and 6. Although a general 
solution to Eq. (3) in the power-series-form can be found, 
two special but very 7 useful cases are considered here, 
which are a pure scaling case (a? 0,1 and P ~0) and a 
pure translation case (a = J and P *0). Not only the 
special cases lead to simpler solutions, but also they are 
more easily implemented in the real measurements. 

The Scaline-Based Method 

In the pure scaling case (a *0,1 and P = <9 ), the 
intrinsic functional equation is 

S(x) = e( X ) - a~‘e( ax ) , ( 5 ) 

where S( x ) is a known function 

S( x ) = m( x ) - a ] m( ax ) . ( 6 ) 

Assume that the functions 8(x) and £( x ) can be 
expanded as a power series in a given domain 
D — [ x } , A ' 2 /, that is, 

N N 

S( x)=^ d n x" and £(x)= ^ e„ x" . ( 7 ) 

n=0 n=0 

For a given set of measurement data points of S( x ) , the 
coefficients d n can be obtained using the least-squares 


method (see Appendix). Substituting Eq. (7) into Eq, (5), 
one can determine all the coefficients e n except for n = I , 

e n = d lt /(l-a n ~ } ). (n=0,2,3-N) (8) 

For n =/, to determine the remaining unknown coefficient 
e l , the measurement value m( x ) is re-written as 

N 

m( x) = ^e„x n +( I + e l )x. (9) 

n=0 

n*i 

Integrating Eq. (9) over a given domain 
D s = [X] S * X 2 X ] ^ & = l x i > x 2 I * onc obtains an 

expression for the coefficient e ] 

2 r* A 

^ = T , I [m(x)-ye n x n Jdx-I, (10) 

*2s ~ X Js “ 

n±l 

where a condition x 2 h * .v j v most hold. We often assign 
/ x !s , x 2s ] = [ x / , / if x ] * x\ . Thus, all the 

coefficients e n are obtained and the bias error distribution 
£( x ) is, in principle, determined. Naturally, the above 
method is called as the scaling-based method for 
recovering the bias error. 

An underlying assumption for the scaling-based 
method is that the functions S( x ) and £(x) can be 
sufficiently expressed as a power series or a polynomial in 
D = [xj ,x 2 ] . To illustrate this method, wc consider a 
hypothetical bias error distribution of a device in 
D = [I.I0] , 

£( x ) = 0.3 — O.Olyfx + 0.03 x -4X10- 4 x 2 -10- 4 X 3 
-6x]0~ 6 -V 5 +0.5 exp( -0.3x 2 ) 

As shown in Fig. 1, the bias error distribution computed by 
using the scaling-based method for a -2 and N = 13 is in 
excellent agreement with one given by Eq. (11). The 
difference between the given and computed bias errors is 
also plotted in Fig. 1 . 

In an ideal case, it seems likely to recover the bias 
error distribution even without using any standard value. 
However, one may notice that both 5( x ) and in Eq. 

(7) are expanded as a function of the true variable x which 
is not exactly known a priori. It will be pointed out later 
that the approximate scaling-based method for actual 
measurements still needs two standard values to recover 
the bias error distribution in the domain D = [x /t x 2 ] . 

The Translation-Based Method 

In the pure translation case ( a-J and P*0 ), 
instead of using a power series, a Fourier series is 
employed because of the shift property of a complex 
exponential function. Since a Fourier series has the 
periodic property, the variable x and the translation 


2 



i 


constant (3 in Eqs. (3) and (4) are replaced by the angular 
variables 0 and G 0 „ respectively. In a domain 
D e = [0,2k] , the intrinsic functional equation for the 
bias error distribution €{6 ) is 

S(0) = e(0)-e(0 + 0 o ), (12) 

where 5(0 ) is a known function for a given 6 0 

S( 6 ) = m( 0 ) - m( 6 + 0 0 ) + 0 o . ( 13 ) 

Assume that the functions 5(6 ) and £( 6 ) in 

D e = [0,2k I can be expanded as a Fourier series 
N 

5(0 ) = ^ [a[ S} cos( nO ) + b\f } sin( nO )] , 

*>=/ 

£(0 ) = a ( 0 e} / 2 + casf j + b[ c) sin( nO ) 7.(14) 

/?=/ 

The coefficients ; and in # ) can be 

determined using the least-squares method for a given set 
of data points (see Appendix). For a relation 

between ( a\f \b f „ €> ) and ( a f „* > , b' n * } ) (n = 1,2, -N) is 
derived by substituting Eq. (14) into Eq. (12), that is, 

(a'/Kb'f'Y =G- , (a , *',h<*') T ,(n = 1.2,3---N) (15) 
where 

( / — cos( n 0 0 ) - sin( n 0 0 ) y 

(j = . 

sin( n0 o ) 1 - cos( n 0 O ) 

Because the determinant det( G ) = 2/ / - cos( n0 0 )] 
should not be zero for the existence of a unique solution, 
the necessary conditions for a unique solution are 
nO a * kn + n / 2 ( k = 0,1,2,-" ) and n * 0 . For n = 0 , 


similar to the pure scaling case, the coefficient a ( 0 c) can be 
determined by the following integral over = [0,2k ] 

/ m( 0 ) 

J () 

N ■ ( 16 ) 

- [ o ( n e } cos( nO ) + b[, c ■ sin( nO ) ]jdO - 2k 

h=i 

Therefore, all the coefficients (a { n c) , b[ e) ) are known and 
the bias error distribution £( 0 ) is readily calculated. We 
refer to this method as the translation-based method. Note 
that the constant term in the Fourier scries for 5(0 ) in Eq. 
(14) is zero. This is a restrictive assumption for the 
translation-based method. Another implicit assumption 
embedded in Eq. (14) is that the bias error distribution 
must be periodic. 

As an example, consider a hypothetical bias error 
£(0 ) defined in D 0 =[0,2k]. 

e( 0 ) = 0.08 + 0.02 cos( 0 )+ 0. 03 cos( 20 ) 

r- - (17) 

- 0.06 sin( 20 )+0.04 t] 0 sin( 40 ) 


Figure 2 shows a comparison between the given bias error 
distribution (17) and the distribution recovered by using 
the translation-based method (0 0 =n/20 and N=6). 
Simulations indicate that the translation-based method is 
good for recovering the bias error distribution whose 
behavior is dominated by trigonometric functions. The 
Fourier series solution is particularly useful for angular 
measurements in the domain D 0 = [0, 2k ] . Compared to 
the scaling-based method, the translation-based method 
has limited applications because of the implicit assumption 
of periodicity 0(0 ) = 0( 2k ). On the other hand, when the 

periodicity condition 0(0 ) = 0( 2k ) is imbedded in data, 
the translation-based method does not explicitly require 
any standard value to recover the bias error distribution. In 
contrast, the scaling-based method typically needs two 
standard values in the real measurements (see Section 3). 

3. The Approximate Scaling-Based Method 

In Section 2, we describe the scaling-based method for 
recovering the bias error distribution as a solution of the 
functional equation. However, the method in the ideal 
condition cannot be directly applied to the real 
measurements because the true variable x in Eq. (7) is not 
known a priori. Therefore, an approximate approach is 
developed, in which the true variable x is replaced by a 
known approximate reference value .v^. The 

approximate reference value x appr is often a measurement 

value by another device with a comparable accuracy. Due 
to the substitution of x appr for a, the scaling-based method 

gives an approximate bias error distribution €(x) appr that 
deviates from the true bias error distribution £(x). The 
deviation of £(x) appr from £(x) can be reasonably 
modeled by a linear function of x uppr . Thus, this 

approximate method requires two standard values to 
determine the unknown coefficients in the linear function 
of deviation for recovering the whole bias error 
distribution. 

Replacing the true variable x in Eq. (7) by a known 
approximate reference variable a- £I/>/ w . , one obtains an 

alternative formulation 

N 

^ ^ n ( ^ a ppr ) 
n=0 

N 

x appr f, (18) 

n=0 

where the function 5(x) appr is an approximate form of 
5( x ) in Eq. (6), i.e., 

S( X = m( X ) - a~’m( a x appr ). (19) 
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The coefficients J n can be obtained by using the least- 

squares method and the coefficients e„ are given by Eqs. 
(8) and (10) in Section 2. Therefore, the approximate bias 
error distribution £( x ) appr is determined. 

Now, we estimate the difference between £(x) and 
£(x) lippr . Assuming x appr = x+£'(x), one knows 


-t hippr ^ ^ n ( ^ appr ) 
n=0 


( 20 ) 


|T 

= £(x)+^n e„ £’( x ) x n ~' + o( £ ’ ) 


n=I 


where C( x ) is the bias error of the approximate reference 
value jr and o(t ) is a higher-order small term of 

e'(x ). Based on the triangle inequality and the Cauchy- 
Schwarz inequality, furthermore, one obtains the following 
estimate for the difference between £( x ) and £( x ) dppr in 

the domain D = [x f ,x 2 } 

W/WeW 


(23) 


As an example, consider a bias error distribution 
e( x ) in D = [ 1, 10 J 

£( x ) = 0.3 - 0.08 Jx -8xl0~ 3 x 2 

+ 8xl0~ 5 x 4 -exp( - 0.3x 2 ) 

The approximate reference value is given by 
Xappr = x + £'( X ) , where 

£’( x ) = A[0.2-5xl0~ 3 4x +0.02x 2 - Jx/0'V 
-0.5 exp(-0.2x) + 0.5 sec h( 0.18.x)] 

The constant A in Eq. (24) is used to adjust the magnitude 
of £'( x ) in computational simulations. Figure 3 shows 
the given bias distribution (23) and the computed 
distributions for a = 2 and N = 13 when the relative 
magnitudes of £ f (x) are II e'W / II £ II = 0.29 , /./7, /.75, 
and 2.95. As indicated in Fig. 3, the computed distribution 
is able to describe the behavior of the bias error 
distribution £( x ) even when the accuracy of the 
approximate reference value is considerably worse than 
that of the tested device. Figure 4 shows the difference 
between the given and computed bias error distributions 


iy 

< (\\£'\\/\\£ 


n=l 


( ..2n-l s .2n-l ) 

1 | x 2 



1/2 

+ C 

(21) 

1 1 £( X J - £( X \ . amp 1 1 / 1 1 £ ( X ) 1 1 

IlfYON/llfTON- The 

as 

linear 

a function of 
relation between 

i 2 ’-‘ j 



\\£(x)-£(x) camp \\/\\£(x)\\ 

and 

W£’(X)W/U€(X)\\ 


where c is a positive constant and 11*11 is the I^-norm showm in Fig. 4 is consistent with the theoretical estimate 


defined as II / 11= [ f f 2 ( x)dx] U2 . Eq. (21) indicates 

that the upper bound of the norm We — £ appr II is simply 
proportional to II £*11 / II £ll , but it is related to the size of 
the domain D = [x Jt x 2 ] in a non-linear fashion. 
Computational simulations indicate that £( x ) appr is often 
a shifted, rotated and sheared representation of £(x) 
although it describes the general behavior of £( x). 

The difference between the true and approximate bias 
error distributions £(x) and £(x) appr can be reasonably 

modeled by a linear function of x appr , that is, 

e( x )~ £( x ) uppr —( C 0 + Cj x appr ) . (22) 

The linear term serves as a correction for the approximate 
solution £(x) appr . Clearly, two standard values are 

required to determine the constants C 0 and C y . When 
two true values are known at two specific points, the 
constants C 0 and C f can be determined. The linear 
model (22) is able to give a reasonable estimate for the 
bias error dislribution. The approximate reference value 
x app/ . is provided by another independent device with the 

accuracy comparable to the tested device. 


( 21 ). 

The selection of the order of the polynomial N in Eq. 
(18) significantly affects the accuracy in calculation of the 
bias error distribution. A polynomial having a low order 
may lead to a poor fit to the true values, while a 
polynomial having an excessively high order may produce 
a large variance due to the over-fitting problem. A 
question is whether there exists an optimal order of the 
polynomial to achieve the highest accuracy of calculation. 
An answer to this problem is not available in a strictly 
mathematical sense. Nevertheless, based on computational 
simulations for various bias error distributions, an 
empirical rule is proposed here for determining an 
appropriate order of the polynomial. We denote £( x f N ) 
as the bias error distribution calculated using a AM-order 
polynomial and define the norm II £( x,N + 1 )-£( x,N )ll 
as the distance between the bias error distributions 
calculated using the (N+f)th- order and iW?-order 
polynomials. Computational simulations show that the 
distance 1 1 £( x t N + I )-£( x,N ) 1 1 usually becomes small 
in a certain range of the order N when the polynomial 
correctly describes the true bias error distribution. This 
phenomenon can be clearly seen when 
We( x,N + J )-£( x,N )\\ is plotted as a function of the 
order N. Figure 5 shows a typical semi-logarithmic plot of 
i\£( x,N +/)-£( x>N ) II as a function of the order N for 
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the given bias error distribution (23), indicating a 
characteristic valley in a range of N = 10-15. In this case, 
a polynomial having the order in N = 10-15 can provide a 
reasonable estimate for the bias error distribution. 

4. Effects of Random Perturbations 

In actual calibration experiments, a measurement 
m(x) has the random error in addition to the bias error. 

The random error affects the solution of the intrinsic 
functional equation for recovering the bias error 
distribution. To simulate effects of the random 
perturbations, we consider a perturbed measurement 
m( x )[ 1 + p( x )] , where p( x) is a random perturbation 

with the normal distribution. Monte Carlo simulations can 
give the probability density distribution of the relative 
difference \\ €( x )-£( x )\\ / \\e( x )\\ that is a random 
variable, where €(x) is the perturbed bias error 
distribution calculated by using the scaling-based method. 
For p(x) having the standard deviation of 0.005, the 
scaling-based method with a = 2 and N= 13 is used to 
recover the bias error distribution (11) for 500 samples. 
Figure 6 shows the probability density distribution of the 
relative difference \\ £(x) — E(x)\\/\\£(x)\\. The 
expectation value of 1! £( x )- £( x )\\ / \\ €( x )\\ is plotted 
in Fig. 7 as a function of the standard deviation of the 
random perturbation p( x ) . 

The scaling constant a may be also changed by a 
small random perturbation in certain measurements. The 
randomly perturbed scaling constant is expressed as 
a( I + p a ) , where p a obeys the normal distribution. 
When p a has the standard deviation of 0.005, the 
disturbed bias error distribution £(x) is calculated by 
using the scaling-based method ( a ~ 2 and N = 13) for 
500 samples for the given bias error distribution (11). 
Figure 8 shows the resulting probability density 
distribution of II £( x)-£( .v)ll / II £(x )\\ , indicating a 
roughly uniform distribution. 

5. Methods for Simultaneously Estimating Bias Errors 
of Two Devices 

The Scalinz-Based Method 

In this section, we describe the scaling-based method 
for simultaneously determining the bias error distributions 
of two devices. Consider devices A and B that have the 
measurements 

m A (x) = x + £ a (x) and m B ( x) = x + e B ( x) . (25) 

The intrinsic functional equations for the bias error 
distributions £ A (x) and £ B (x) arc 

Sj(x) = £ a (x)-£ r (x) 

" , (26) 

5 2 (x) = £ a (x)~ a ! £ b (ax) 


where 8 j( x ) = m A ( x ) - m B ( x) and 

S 2 ( x ) = m A ( x ) - a~ I m B ( ax ) arc known functions 
obtained from measurements. Assume that the functions 
8j( x ) , S 2 ( x) , £ a ( x ) and £ B ( x) can be expanded as a 

power scries in a given domain D ~ / x } , x 2 / , that is, 

N N 

S } ( x) = ^<I ( „ n x n and S 2 ( x) = ^d[ 2 V 1 , 

n-0 n-0 

N N 

£,\( x) = y ^e ( „ A) x n and £„(x) = . (27) 

n~0 n-0 

For a given set of data points of 5 { ( x) and S 2 (x), the 

coefficients d ( J * and d[ 2) can be obtained using the 
least-squares method. Except for n -I, the coefficients 
e\ } A) and e { n B } can be determined by 


d i 2 ) -a"‘id (n 
el A >= " , " and 


d t2 >-d (n 

( B) _ u n u n 


]-a n ~ } ' i-a n ~ ] 

(n =0,2,3 — N) (28) 

For n =/, since the functional equations (26) are reduced 
to one independent equation, there is only a single 

e ( } B * for two unknowns 


i 


equation cl y = d) ' = e 
e { j A) and e \ B} . To determine the remaining unknowm 


coefficients e) and e) , we use an integral expression 

N 

e \ A} = J- 2 I dx - J , (29) 

x 2s ~ *ls Jx » 

n*f 


and the relation e (B) -e ( , A) -d\ n . The given domain 
D A = [ x lst x 2s ] ( xj y & x 2s ) for integration is a sub- 
domain of the measurement domain D = [ x h x 2 ] . 

As discussed previously, since the true variable v is 
not known a priori in practical applications, approximation 
should be used in which the true variable x is replaced by 
the known approximate variable m A (x). The 
approximate formulations are 

N 

5,(x) appr =^d ( n n lm A (x)r, 

n-0 

yv 

S 2 (x) appr =X d ( n 2, f m A ( x)]" , 

n=0 

N 

Ca( x U w m ,\(x)]\ 

n-0 

£/t( x ) appr = ^ e\, B 1 1 m A ( x )] ” . (30) 

n=0 
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where 8 } (x) appr =m A (x)-m B (x) and 

S 2 (x) appr ^m A (x)-a^m it (am A (x)) are known 

functions. The approximate bias errors £ A (x) appr and 

e B ( x ) uppr can he determined by using Eqs. (28) and (29). 

The approximations C A (x) uppr and £ B ( x) appr describe 

the general behavior of the bias error distributions, but 
deviate from the true ones. The linear model for correcting 

£a( x Uw and € *(- x kw is s ivenb y 

£ A ( = ^ A ( ^ )uppr ~~ f ^ AO + A I ^ A^ ' ^ 

£ B (x) = £ B (x ■) appr -[C H0 +C BJ m A (jc)} . (31) 

When two true values are known at two specific points, the 
constants G AB * ^ai » ^bo * and G B j can be determined. 

Simulations indicate that Eq. (31) is able to describe a 
class of the bias error distributions £ A ( x) and £ B (x) that 
arc reasonably represented by a power series. As an 
example, we consider the bias error distributions £ A ( x) 

and £ B (x) in D = [IJ()J 

e A (x ) = 0.3 -4x 10 ~ 3 Jx - 2 x W 4 x 2 - 8 x 10~- x 4 , 

e B ( x ) = 0./-5x 10~* -Jx + 2 x I0~ 2 x 2 - 3x 10~ 4 x J . (32) 
Figure 9 shows the given and calculated bias error 
distributions £ A (x) and e B (x) for a = 2 and N = 10. 
The empirical rule for selecting the order of the 
polynomials in Section 3 is also applicable to the two- 
device case. 

The Translation-Based Method 

In a domain D 0 = [0 } 2/r ] , the intrinsic functional 

equations for the bias errors s A ( 9 ) and e B ( 9 ) are 

S,(0)-e A (0)-e B <0) 33 

S 2 (0) = £ a ( 0)~ £g(0+0 o )’ 

where 9 0 is a constant translation in radian, 
S } (9 )- m A ( 9 )-tn B ( 9 ) , and 

8 2 (0 ) = ni A (0 )-m B (0 + 0 o ) + 6 0 are known functions. 
Assume that in D o =[0,2n] the functions 8 } (0), 
S 2 (0)> £ a (9) and £ B (9) can be expanded as 

N 

5 l (8) = a l 0 n /2 + ^la < „ n cos( n 6 ) +b[ ' 1 sin( n8)J, 

j V 

S 2 ( 0 ) = a { 0 2 ] / 2 + l a ( „ 2 } cos( n 9 ) + h[ 2i sin( n9 )} , 

fi—i 

N 

£ A ( 6 ) = a ( 0 A V 2 + ^ / a ( fi A ) cos( ti9 )+h ( „ A } sin( n9 )] , 

n=! 

N 

e B ( e ) = a ( 0 B ■ / 2 + £ I a \; 8 ' cos( n 9 ) 4 b ( n B ; sin{ n9)J. 

i i=i 

(34) 


The coefficients (a { n 1} , b ( fl 1 } , a[ 2) , b ( ti 2) ) in 8 { (9) and 
5i(9 ) can be determined using the least-squares method. 
For n*0 } a relation between {a ( n A) , b ( n A) , ct (B) , b ( t Br ) 
and ( a \,' 1 , b ( „’ 1 , a[ 2 J , b [ 2 ' ) is 

(<’ «»' try =c-'(«"> hir air birY . 

( n = 1,2,3 ■■■ N ) (35) 

where 

f l 0 -1 0 y 

0 10 -1 
G — 

1 0 -cos(n8 0 ) -sin(n9 0 ) 

1 0 1 sin(n0 o ) -cos(n9 0 ) ^ 

Since the determinant of G is det( G )-2f 1 - eos( ti0 o )/ , 
a necessary condition for the existence of a unique solution 
is n6 0 * kn + 7T / 2 ( k =0,/,2,-” ). For u-0, only one 


equation a ( Q ,} = a{ 2) - a ( 0 A) -a t 0 B> is available for two 
(B) 


unknowns a ( 0 A) 


and a f 0 l1) - Using a similar method to the 


pure scaling case, the coefficient a ( 0 A) can be determined 
by an integral over D e = [0 y 27t ] 

2k 


a u Ai =tt 'I (m A (0) 

0 

N 

- y la' n A 1 cost nOj+b ( „ A) sin( n 0 ) ]}d0 - 2n 


(36) 


H = 1 

Therefore, (he coefficient a l 0 R> - a' (J A> - a { Q n is readily 
known. As an example, we consider the following bias 
errors e A ( 8 ) and e B ( 6 ) in D g = [ 0,2n ] 

£ a (8) = 0.03 + 5x 10- 4 8 + 0.02 cas( 8) + 3x 10 ~ * cos( 28 ) 
-6 x /(T ? sin( 8 ) + 4xl0~ 4 sin( 38 ) 
e H (0 ) = 0.05 + 3x 10 4 8 + 0.05 cos( 8)- 6x10^ cos( 38 ) 
+ I0- 4 sin( 8 ) + 3x10-' sin( 48 ) 


(37) 

Figure 10 shows a comparison between the given bias 
errors and those computed by using the translation-based 
method ( 0 0 = n/ 35 and N -6). 


6 . Applications 

In principle, the aforementioned methods are 
applicable to any device. Here, to illustrate applications of 
these methods, we present several typical examples in 
voltage measurements, angular measurements, and optical 
measurements. 

Volta ge Calibrations for A/D Converter 

A voltage divider was constructed using stable 
resistors to obtain a nominal voltage ratio of 0.5 (the 
scaling constant 0.5 in the scaling-based method). The 
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input to the voltage divider was connected to a stable 
voltage source and the output was connected to one 
channel of a 24 bit multiplexing A/D converter (Lawson 
Model 201). Another channel was connected directly to 
the voltage source. The voltage source was varied from - 
5v to +5v (the full range of the A/D converter). The 100 
readings on both channels were recorded. No accuracy 
was assumed for the voltage source or the voltage divider 
ratio. Calculation for the bias error distribution was 
performed by using the approximate scaling-based method. 
Two reference readings were taken to provide the standard 
values, one with the input shorted and another with the 
input voltage set to exactly 3 volts (as measured on a 
precision voltmeter (HP 3458A) of the accuracy of 8ppm). 
This is all that is required to find the calibration of the A/D 
converter by using the approximate scaling-based method. 
Eleven additional points were taken using the precision 
voltmeter (HP 3458A) over the full range to verify the 
mathematical solution. The voltage calibration plots in 
Fig. 1 1 shows the verification points and the calibration 
results computed using the scaling-based method for a — 
0.5 and N ~ 15. The calibration curve computed by using 
the scaling-based method is in good agreement with the 
verification data. 

Instead of using a voltage divider for scaling, scaling 
can be achieved based on dial readings of a power supply 
(Wavelek Model 220). An advantage of this approach is 
that a physical f device like the divider is no longer required 
for scaling. Figure 1 1 shows the bias error distribution 
computed by using the scaling-based method in which 
scaling is controlled based on the dial readings. Compared 
with the results given by using the divider, the distribution 
obtained by using the dial readings shows larger local 
variations. This is because the mechanical voltage setting 
device of the power supply caused abrupt voltage 
variations over a certain range of operation. It is expected 
that a more stable device will produce a smoother bias 
error distribution. 

Angle Calibrations for Encoder and Indexing Table 

The two-device translation-based method was used for 
angle calibration on a dividing head with attached encoder. 
The secondary device was an indexing table with a 
resolution of one degree. Both devices have a nominal 
accuracy specification of one arc second. The axis of 
rotation was horizontal so that a precision servo 
accelerometer was used as an indicator of level. The 
procedure involves rotating the dividing head clockwise in 
10° increments, rotating the indexer counter-clockwise in 
10° increments, and reading the level as indicated by the 
precision accelerometer. The second set of data, taken 
with the indexer and accelerometer translated -60 degrees. 
To recover the bias error distributions, the first set of data 
and the second set of -60°-tran slated data were processed 
by using the two-device translation-based method with a 
4th-order Fourier scries. The calibration curves for both 
the encoder and indexing table can be simultaneously 


determined, as shown in Fig. 1 2. The black circular points 
in Fig. 12 are the results of a conventional calibration 
using a device with an accuracy four times belter than the 
dividing head with the encoder. The bias error distribution 
measured by using the precision accelerometer is in good 
agreement with the computed distribution for the encoder. 
It is worthwhile noting that we do not explicitly use any 
standard value to recover the bias error distribution. 
Nevertheless, the periodicity condition 6(0)- 9(271), 
which is automatically satisfied in the Fourier scries 
solution, can be considered as an imposed constraint. The 
two-device translation-based method, originally proposed 
by Finley [31, has been used regularly by Wyle 
Laboratories to calibrate precision angle-measuring 
devices in their facility and at NASA Langley. 

Radial Lens Distortion 

An interesting example is application of the scaling- 
based method to determination of the radial lens distortion. 
In reality, a lens used for imaging is not perfect and the 
imperfect lens may distort an image. Thus, camera 
calibration to determine the camera parameters including 
the lens distortion parameters is crucial for accurate image- 
based measurements [7). The most dominant lens 

distortion is the radial lens distortion that is symmetric 
about the principal-point (close to the geometric center of 
an image). The radial lens distortion is described by a 

simple model Sr = K { r 3 , where Sr is a bias error in the 
radial distance due to the lens distortion, K } is the radial 
distortion parameter, and V is the radial distance from the 
principal-point. The radial distortion parameter in 
addition to other camera parameters can be obtained in 
comprehensive camera calibrations by using analytical 
photogrammetric techniques. For an 8mm Computar TV 
lens used in this test, an optimization camera calibration 

method [7] gives K } — 0.00 J 297 mm' 2 . 

Unlike analytical photogrammetric techniques that use 
a special mathematical model for the lens distortion, the 
scaling-based method determines the radial lens distortion 
under a general theoretical framework of the bias error. 
During tests, a 22inxl7in target plate with 121 retro- 
reflective targets of l/2in diameter was used to provide a 
planar target field. A CCD video camera (Hitachi KP- 
FIU) with an 8mm Computar TV lens, viewing 
perpendicularly the target plate, was used to lake images of 
the plate placed at two different distances from the camera. 
Figure 13 shows typical images of the target plate at tw'o 
different distances from the camera. The digitized image 
has 640x480 pixels and the nominal pixel spacing is 9.9 
pm in both the horizontal and vertical directions. The 
centroids of the targets and the radial distances of the 
targets from the geometric center in these images were 
computed. Two images of the target plate at two different 
distances from the camera naturally provide scaling in the 
radial distance in image plane. The scaling constant 
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a = 0.8262 was obtained by averaging ratios between the 
radial distances of the corresponding targets from the 
geometric center in two images. The bias error 
distribution £( r ) was calculated by using the approximate 
scaling-based method when the order of the polynomial is 
suitably chosen in N = 20-28. By the definition of e( r ), 
one knows Sr( r) = -e( r ) . The condition Sr(0) = 0 
must be satisfied. Thus, only one standard value is 
required to determine the unknown constants in Eq. (22), 
which is given by the optimization camera calibration 
method. Figure 14 shows the radial lens distortion 
distributions obtained by the scaling-based method (N = 
25) and the optimization camera calibration method for an 
8mm Computar TV lens. 

7. Conclusions 

The bias error distribution of a device can be sought as 
a solution of the intrinsic functional equation. Based on 
this idea, the scaling- and translation-based methods have 
been developed to determine the bias error distribution in a 
domain from less accurate measurements. The scaling- 
based method is applicable to a device whose bias error 
distribution can be adequately expressed as a power series 
or a polynomial in a domain. Practical application of the 
sealing-based method typically requires two standard 
values to recover the complete bias error distribution in the 
whole domain. The suitable order of a polynomial for 
accurate recovery of the bias error distribution can be 
selected according to an empirical rule proposed in this 
paper. The translation-based method, which uses a Fourier 
series to describe the bias error distribution, is particularly 
useful for angular measurements because of the internal 
periodicity constraint. The translation-based method does 
not explicitly require any standard value. These methods 
have been extended to simultaneously determine the bias 
error distributions for two devices. To validate and clarify 
the technical aspects of these methods, computational 
simulations have been carried out for various hypothetical 
distributions of the bias error. Laboratory tests have, been 
conducted for calibrating several different devices such as 
A/D converter, angular measurement device and optical 
lens. These methods of estimating the bias error 
distribution are effective for a variety of devices. 
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Appendix: Least-Squares Estimation of the Coefficients 

To determine the coefficients cl n (n = 0,1,2, ••• N ) 
in the power series or polynomial S( x ) in Eq. (7) for a set 
of data points x t (i = 1,2, M ), a system of equations 
for d n ( n =0,1,2, N ) is 

Pd =.<L (AI) 

where d = {cl () d } cl ? • * • d N ) , 

S = (S(x, ) S(x 2 ) ■■■ S(x M )) T , and 

ft 2 2 N \ 

I X, XJ XJ ••• X, 

I 2 2 N 

p_ X 2 X 2 x 2 x 2 

i 2 .2 N 

J X M X M X M *** X M , 

The least-squares solution to (A I ) is 

d = (P T P) 1 P T 5 . (A2) 

Similarly, a system of equations for the coefficients 
( a ( n 5} , b { u S} } in the Fourier scries in Eq. (14) for a set of 
data points 0 f ( i = /, 2, • • • M ) is 

Fa = S , (A3) 

wherea = \r/ / ' mmm 0 N bj '”b N j . 

S-{S(0j )S(0 2 ) )) t , and 

f cos( Gj ) • • ■ cost N6 1 ) sin( 6 } )■■- sm( N6 } ) ^ 

cos( 6 1 ) ♦ • • cos( N 6 7 ) sin( 0 7 ) • • * sin( NG 7 ) 

F - " " ' 

^cos( 0 M ) • ■ • cos( N 0 M ) sin( 0 M )••• sin( N 0 M ) ^ 
The least-squares solution to (A3) is 

a = (F r F)' 1 F T S . (A4) 
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Fig. 1 . The given bias error distribution and computed 


results by using the scaling-based method (a= 2, N = 13). 



Fig. 2. The given bias error distribution and calculated 
results by using the translation-based method 
(0 0 =71/20, N =6 ). 



X 


Fig. 3. The given bias error distribution and calculated 
results by using the approximate scaling-based method 
(a=2,N= 1 3) for four magnitudes of the bias error of the 
approximate reference value. 



Ilf’{ A )ll/llf{ V >n 

Fig. 4. The difference between the given and computed 
bias errors as a function of 1 1 £'( x ) 1 1 / 1 1 S( x ) 1 1 . 
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Fig. 6. The probability density function of 
1 1 €( x ) - £( x ) i i / 1 1 e ( x ) 1 1 for a randoml y perturbed 
measurement. 
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std of perturbation 


Fig. 7. The expectation value of II 2( x ) — €( x )\\ / \\ €{ x )\\ 
as a function of the standard deviation of the random 
perturbation. 



II £( -V ) — .V ^ II / II f{ .V ) II 

Fig. 8. The probability density function of 
II £( x ) - e( x ) 1 1 / 1 1 £( x ) 1 1 for a randomly perturbed 
scaling constant. 



Fig. 9. The given bias error distributions and computed 
results by using the approximate two-device scaling-based 
method ( <7= 2, N - 10) for devices A and B. 
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Fig. 10. The given bias error distributions and computed 
results by using the two-device translation-based method 
( 6 0 = ; i/ 35 , N = 6 ) for devices A and B. 



Voltage (volt) 

Fig. 1 1. The measured bias error distribution of voltage 
and computed results by using the scaling-based method 
(a= 0.5, N = 15) for an A/D converter. 



Angle (radian) 

Fig. 12. The measured bias error distribution of angle and 
computed results by using the two-device translation-based 
method ( 0 0 = -it / 3 , N = 4 ) for an encoder and an 
indexing table. 


10 








Fig. 14. The radial lens distortion distribution computed by 
using the scaling-based method (a= 0.8262, N=25) for of 
an 8mm Compuiar TV lens, along with the camera 
calibration result obtained by using an optimization 
method. 


Fig. 1 3. Images of the target plate at two different 
distances relative to the camera. 
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